Clase 2: Operaciones Espaciales
Analisis espacial con R
Martín Alalú
29/04/2021
1 Operaciones espaciales
Hay una enorma cantidad de operaciones que se le pueden hacer a los objetos geográficos. Calcular el área, el perímetro, el largo, distancias entre objetos, intersecciones, uniones, diferencias, realizar transformaciones y más! Estas manipulaciones son sumamente útiles ya que nos permitirán comprender mejor el fenómeno que estamos analizando.
Para más información de todas las operaciones posibles (o las principales) sugerimos revisar el capítulo 4 del libro Geocomputation with R.
2 Área, perímetro e intersección
Como siempre comenzamos cargando las librerías que vamos a usar. Recuerden que si no las tienen instaladas aún lo pueden hacer con install.packages(“NOMBRE_LIBRERIA”).
De paso cargamos un dataframe geográfico (formato .shp) que tiene los polígonos de las provincias de Argentina.
library(sf)
library(tidyverse)
library(units)
options(scipen = 999)
prov <- read_sf("https://raw.githubusercontent.com/martoalalu/clase-geo-salud/clase-2021/data/ar_provincias.geojson")Veamos qué información tiene.
head(prov)## Simple feature collection with 6 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -72.13788 ymin: -46.0069 xmax: -56.66736 ymax: -24.09314
## Geodetic CRS: WGS 84
## # A tibble: 6 x 2
## provincia geometry
## <chr> <MULTIPOLYGON [°]>
## 1 Buenos Aires (((-62.4578 -38.75891, -62.45403 -38.75514, -62.43847 -38.…
## 2 Catamarca (((-64.82775 -29.56298, -64.9697 -29.58564, -64.90726 -29.…
## 3 Chaco (((-60.9271 -27.99942, -61.71717 -27.99393, -61.71229 -25.…
## 4 Chubut (((-66.36792 -45.04317, -66.37458 -45.05514, -66.36292 -45…
## 5 Ciudad de Buenos … (((-58.44724 -34.69322, -58.46166 -34.70905, -58.52893 -34…
## 6 Córdoba (((-63.39265 -34.99844, -65.10672 -34.99768, -65.10981 -33…Ok, sólo el nombre y el campo de geometría, obvio. Hagamos un mapa!
ggplot(prov) +
geom_sf()
Podemos rellenar cada provincia con un color distinto.
ggplot(prov) +
geom_sf(aes(fill=provincia))
Ok. Ahora sí, empecemos a agregarle información a este dataframe. Calcular el área es muy simple con la función st_area de la librería sf.
st_area(prov) #nos devuelve en unidades raras, pasamos a km2## Units: [m^2]
## [1] 306889486297 101887635888 99886544213 224258973324 209979717
## [6] 164683302391 88973340510 77861018534 75596524598 53192452823
## [11] 142688382924 90990885271 148787352879 29982817814 94604326333
## [16] 202590178533 155322607325 88301731606 75855929854 243981330500
## [21] 132921440646 136871676316 20909206585 22537877676Mmm suenan raros estos números. Pasemos la unidad de medición a kilómetros cuadrados!
set_units(st_area(prov), km^2)## Units: [km^2]
## [1] 306889.4863 101887.6359 99886.5442 224258.9733 209.9797 164683.3024
## [7] 88973.3405 77861.0185 75596.5246 53192.4528 142688.3829 90990.8853
## [13] 148787.3529 29982.8178 94604.3263 202590.1785 155322.6073 88301.7316
## [19] 75855.9299 243981.3305 132921.4406 136871.6763 20909.2066 22537.8777Ahí va mejor. Pero seguimos sin saber a qué provincia corresponde cada registro. Lo pasamos a formato numérico y agregamos como columna nueva!
prov$area <- as.numeric(set_units(st_area(prov), km^2))
prov## Simple feature collection with 24 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -73.57778 ymin: -55.05736 xmax: -53.59184 ymax: -21.77855
## Geodetic CRS: WGS 84
## # A tibble: 24 x 3
## provincia geometry area
## * <chr> <MULTIPOLYGON [°]> <dbl>
## 1 Buenos Aires (((-62.4578 -38.75891, -62.45403 -38.75514, -62.4384… 3.07e5
## 2 Catamarca (((-64.82775 -29.56298, -64.9697 -29.58564, -64.9072… 1.02e5
## 3 Chaco (((-60.9271 -27.99942, -61.71717 -27.99393, -61.7122… 9.99e4
## 4 Chubut (((-66.36792 -45.04317, -66.37458 -45.05514, -66.362… 2.24e5
## 5 Ciudad de Bueno… (((-58.44724 -34.69322, -58.46166 -34.70905, -58.528… 2.10e2
## 6 Córdoba (((-63.39265 -34.99844, -65.10672 -34.99768, -65.109… 1.65e5
## 7 Corrientes (((-57.78307 -30.43265, -57.84768 -30.47401, -57.867… 8.90e4
## 8 Entre Ríos (((-58.5247 -33.46152, -58.53116 -33.49211, -58.5254… 7.79e4
## 9 Formosa (((-58.24602 -26.73134, -58.25027 -26.75323, -58.261… 7.56e4
## 10 Jujuy (((-64.83212 -24.47771, -64.86935 -24.5232, -64.8812… 5.32e4
## # … with 14 more rowsBien! Ahora que tenemos una columna con el área podemos hacer un mapa coloreando cada provincia según los km2.
ggplot(prov) +
geom_sf(aes(fill=area))
Mmm está mal, pero no tan mal. Lo que ya sabíamos, la Provincia de Buenos Aires es la más grande del país, pero también se ve cierto patrón, a medida que pasamos de Norte a Sur el gradiente se ve más claro, o sea que las provincias suelen ser más grandes…
Hagamos un gráfico de barras para ver este patrón.
ggplot(prov) +
geom_bar(aes(x= reorder(provincia, area), weight=area, fill=area)) +
coord_flip() 
Claro, al pasar del mapa al gráfico de barras perdemos la referencia geográfica. Agreguemos una columna con la región de la provincia para ver si las de la Patagonia suelen ser más grandes que las del Norte.
prov <- prov %>%
mutate(region =
case_when(
provincia %in% (c("Buenos Aires", "Córdoba", "La Pampa", "Entre Ríos", "Santa Fe", "Ciudad de Buenos Aires")) ~ "Pampeana",
provincia %in% (c("Mendoza", "San Luis", "San Juan")) ~ "Cuyo",
provincia %in% (c("La Rioja", "Catamarca", "Jujuy", "Salta", "Tucumán", "Santiago del Estero")) ~ "Noroeste",
provincia %in% (c("Corrientes", "Misiones", "Formosa", "Chaco")) ~ "Noreste",
TRUE ~ "Patagonia"
)
)
ggplot(prov) +
geom_bar(aes(x= reorder(provincia, area), weight=area, fill=region)) +
coord_flip() 
Ahí va mejor! Hay cierto patrón, 3 de las 4 provincias más grande son de la Patagonia y 5 de las 6 más chicas del Norte! Esto nos sirve para entender que al hacer análisis espacial a veces salir del mapa puede ser sumamente útil para ver patrones.
Una más y listo. Agreguemos una línea vertical con la media del área de las provincias para ver cuáles están por encima y cuáles por debajo.
ggplot(prov) +
geom_bar(aes(x= reorder(provincia, area), weight=area, fill=region)) +
geom_hline(yintercept=mean(prov$area), linetype="dashed",
color = "black", size=0.5) +
coord_flip() 
¿Qué conclusiones sacan ahora?
Bueno, sigamos con el perímetro. La función es st_length.
st_length(prov)## Units: [m]
## [1] 4502992.56 2045513.65 1787796.08 3069110.91 61337.54 1887257.68
## [7] 1523652.60 1371928.96 1823365.32 1356250.58 1801249.30 1592326.16
## [13] 2012790.78 1218772.88 1767529.09 2570198.77 2869031.24 1716796.00
## [19] 1364124.16 2773527.11 2143460.03 1701827.16 1449042.72 806231.04Vemos que la unidad de medida son metros, lo pasamos a km con solo dividirlo por 1000 y agregamos a columna nueva.
prov$perimetro <- as.numeric(st_length(prov) / 1000)
prov## Simple feature collection with 24 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -73.57778 ymin: -55.05736 xmax: -53.59184 ymax: -21.77855
## Geodetic CRS: WGS 84
## # A tibble: 24 x 5
## provincia geometry area region perimetro
## * <chr> <MULTIPOLYGON [°]> <dbl> <chr> <dbl>
## 1 Buenos Aires (((-62.4578 -38.75891, -62.45403 -38.7… 3.07e5 Pampe… 4503.
## 2 Catamarca (((-64.82775 -29.56298, -64.9697 -29.5… 1.02e5 Noroe… 2046.
## 3 Chaco (((-60.9271 -27.99942, -61.71717 -27.9… 9.99e4 Nores… 1788.
## 4 Chubut (((-66.36792 -45.04317, -66.37458 -45.… 2.24e5 Patag… 3069.
## 5 Ciudad de Bu… (((-58.44724 -34.69322, -58.46166 -34.… 2.10e2 Pampe… 61.3
## 6 Córdoba (((-63.39265 -34.99844, -65.10672 -34.… 1.65e5 Pampe… 1887.
## 7 Corrientes (((-57.78307 -30.43265, -57.84768 -30.… 8.90e4 Nores… 1524.
## 8 Entre Ríos (((-58.5247 -33.46152, -58.53116 -33.4… 7.79e4 Pampe… 1372.
## 9 Formosa (((-58.24602 -26.73134, -58.25027 -26.… 7.56e4 Nores… 1823.
## 10 Jujuy (((-64.83212 -24.47771, -64.86935 -24.… 5.32e4 Noroe… 1356.
## # … with 14 more rowsExcelente! Ahora sí, al mapa.
ggplot(prov) +
geom_sf(aes(fill=perimetro))
Y al gráfico de barras.
ggplot(prov) +
geom_bar(aes(x= reorder(provincia, perimetro), weight=perimetro, fill=region)) +
geom_hline(yintercept=mean(prov$perimetro), linetype="dashed",
color = "black", size=0.5) +
coord_flip() 
3 Intersecciones
Otra de las operaciones espaciales más útiles es la de calcular intersecciones entre 2 objetos espaciales. Por ejemplo ver qué puntos (escuelas) caen dentro de un polígono (ej. barrio). Muchas veces la información disponible es de un universo mucho mayor a nuestro fenómeno de estudio y necesitamos hacer un recorte, ahí es dónde st_intersects nos puede ayudar.
Asimismo por ahora estuvimos manejando un tipo de dato geográfico (polígonos) y haciendo un estilo de mapas (coropléticos), probemos usando datos que sean puntos y hagamos otros tipos de mapas.
Vamos a usar un dataset disponibilizado por el Ministerio de Salud del Gobierno Nacional que contiene la ubicación exacta todos los efectores de salud del país.
#Cargamos los datos
salud <- read.csv("https://github.com/martoalalu/clase-geo-salud/raw/master/data/refes-hospitales.csv", fileEncoding = 'UTF-8')
#Vemos la clase de objeto que es
class(salud)## [1] "data.frame"head(salud)## CODIGO GLN NOMBRE
## 1 50380842150137 NA PUESTO EL TORO (17)
## 2 50380632150075 NA PUESTO BARRO NEGRO (9) (EX LOTE ROSARIO DE RIO GRANDE)
## 3 50260212329179 NA KIN SPORT
## 4 51260352329290 NA LABORATORIO ANALISIS CLINICOS GEROSA
## 5 51260352329289 NA LABORATORIO ANALISIS CLINICOS
## 6 51260352329293 NA LABORATORIO ANATOMIA PATOLOGICO (SALVAGNO)
## CODIGO_TIPOLOGIA
## 1 ESSIDT
## 2 ESSIDT
## 3 ESSIDT
## 4 ESSID
## 5 ESSID
## 6 ESSID
## TIPOLOGIA
## 1 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 2 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 3 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 4 Establecimiento de salud sin internación de diagnóstico
## 5 Establecimiento de salud sin internación de diagnóstico
## 6 Establecimiento de salud sin internación de diagnóstico
## CODIGO_CATEGORIA_TIPOLOGIA
## 1 ESSIDT
## 2 ESSIDT
## 3 ESSDIT
## 4 ESSID
## 5 ESSID
## 6 ESSID
## CATEGORIA_TIPOLOGIA
## 1 Sin atención médica en forma periódica (menor a 3 veces por semana)
## 2 Con atención médica general por lo menos 3 días de la semana
## 3 Con atención médica diaria y con especialidades y/o otras profesiones
## 4 Laboratorio de Análisis Clínicos
## 5 Laboratorio de Análisis Clínicos
## 6 Laboratorio de Anatomía patológica
## CODIGO_DEPENDENCIA DEPENDENCIA ORIGEN_FINANCIAMIENTO CODIGO_PROVINCIA
## 1 21 Provincial Público 38
## 2 21 Provincial Público 38
## 3 23 Privado Privado 26
## 4 23 Privado Privado 26
## 5 23 Privado Privado 26
## 6 23 Privado Privado 26
## PROVINCIA CODIGO_DEPARTAMENTO DEPARTAMENTO CODIGO_LOCALIDAD
## 1 Jujuy 84 Susques 38084030000
## 2 Jujuy 63 San Pedro 38063090000
## 3 Chubut 21 Escalante 26021030018
## 4 Chubut 35 Futaleufú 26035030000
## 5 Chubut 35 Futaleufú 26035030000
## 6 Chubut 35 Futaleufú 26035030000
## LOCALIDAD CP DOMICILIO TE1
## 1 EL TORO 4411 AVENIDA LAVALLE ESQUINA 24 DE OCTUBRE
## 2 LA MENDIETA 4522 CARLOS GARDEL S/N° - LOTE BARRO NEGRO
## 3 COMODORO RIVADAVIA 3240 Avenida Rivadavia 633 0297-4475310
## 4 ESQUEL 9200 25 de Mayo 498 2945452050
## 5 ESQUEL 9200 Ameghino 1033 02945-452155
## 6 ESQUEL 9200 Belgrano 347
## TE2 TE3 TE4 LATITUD LONGITUD ESTADO_GEO
## 1 -23.20853 -66.82709 SINCONFIRMAR
## 2 -24.30814 -64.93175 SINCONFIRMAR
## 3 -45.86213 -67.48355 CONFIRMADO
## 4 294515503809 -42.91505 -71.31971 CONFIRMADO
## 5 -42.91577 -71.31803 CONFIRMADO
## 6 -42.91788 -71.31992 SINCONFIRMARsummary(salud)## CODIGO GLN NOMBRE
## Min. :10020012215200 Min. :7798100490010 Length:28366
## 1st Qu.:50063572202800 1st Qu.:9990454430010 Class :character
## Median :50500072360600 Median :9991301000120 Mode :character
## Mean :43797418878000 Mean :9868266292960
## 3rd Qu.:51065682305200 3rd Qu.:9992051100060
## Max. :53940142795200 Max. :9992574600000
## NA's :28241
## CODIGO_TIPOLOGIA TIPOLOGIA CODIGO_CATEGORIA_TIPOLOGIA
## Length:28366 Length:28366 Length:28366
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## CATEGORIA_TIPOLOGIA CODIGO_DEPENDENCIA DEPENDENCIA
## Length:28366 Min. :20.00 Length:28366
## Class :character 1st Qu.:22.00 Class :character
## Mode :character Median :23.00 Mode :character
## Mean :22.52
## 3rd Qu.:23.00
## Max. :32.00
##
## ORIGEN_FINANCIAMIENTO CODIGO_PROVINCIA PROVINCIA CODIGO_DEPARTAMENTO
## Length:28366 Min. : 2.00 Length:28366 Min. : 1.0
## Class :character 1st Qu.: 6.00 Class :character 1st Qu.: 28.0
## Mode :character Median :30.00 Mode :character Median : 70.0
## Mean :38.81 Mean :157.5
## 3rd Qu.:66.00 3rd Qu.:140.0
## Max. :94.00 Max. :882.0
##
## DEPARTAMENTO CODIGO_LOCALIDAD LOCALIDAD CP
## Length:28366 Min. : 6028010 Length:28366 Length:28366
## Class :character 1st Qu.: 6588100000 Class :character Class :character
## Mode :character Median :30035045000 Mode :character Mode :character
## Mean :37490327544
## 3rd Qu.:66028050000
## Max. :94014020000
##
## DOMICILIO TE1 TE2 TE3
## Length:28366 Length:28366 Length:28366 Length:28366
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## TE4 LATITUD LONGITUD ESTADO_GEO
## Length:28366 Min. :-54.83 Min. :-72.89 Length:28366
## Class :character 1st Qu.:-34.70 1st Qu.:-65.35 Class :character
## Mode :character Median :-32.91 Median :-61.93 Mode :character
## Mean :-32.52 Mean :-62.47
## 3rd Qu.:-28.47 3rd Qu.:-58.67
## Max. :-21.94 Max. :-53.65
## NA's :2240 NA's :2240Lo cargamos como un dataframe normal no-geográfico pero al explorarlo vemos que tiene 2 columnas que sí precisan su ubicación exacta en la tierra, LATITUD y LONGITUD. Tenemos que hacerle entender a R que se trata de un dataframe geográfico, es fácil!
Primero filtramos y sólo nos quedamos con las observaciones que no son nulas ya que no podemos ubicar en un mapa objetos con latitud y longitud nulas. Luego con st_as_sf convertimos el dataframe a un objeto geográfico indicandole dónde están los datos de las coordenadas y le indicamos el sistema de referencia de coordenadas.
salud <- salud %>%
filter(!is.na(LONGITUD), !is.na(LATITUD)) %>%
st_as_sf(coords = c("LONGITUD", "LATITUD"), crs = 4326)
class(salud)## [1] "sf" "data.frame"head(salud)## Simple feature collection with 6 features and 23 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -71.31992 ymin: -45.86213 xmax: -64.93175 ymax: -23.20853
## Geodetic CRS: WGS 84
## CODIGO GLN NOMBRE
## 1 50380842150137 NA PUESTO EL TORO (17)
## 2 50380632150075 NA PUESTO BARRO NEGRO (9) (EX LOTE ROSARIO DE RIO GRANDE)
## 3 50260212329179 NA KIN SPORT
## 4 51260352329290 NA LABORATORIO ANALISIS CLINICOS GEROSA
## 5 51260352329289 NA LABORATORIO ANALISIS CLINICOS
## 6 51260352329293 NA LABORATORIO ANATOMIA PATOLOGICO (SALVAGNO)
## CODIGO_TIPOLOGIA
## 1 ESSIDT
## 2 ESSIDT
## 3 ESSIDT
## 4 ESSID
## 5 ESSID
## 6 ESSID
## TIPOLOGIA
## 1 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 2 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 3 Establecimiento de salud sin internación de diagnóstico y tratamiento
## 4 Establecimiento de salud sin internación de diagnóstico
## 5 Establecimiento de salud sin internación de diagnóstico
## 6 Establecimiento de salud sin internación de diagnóstico
## CODIGO_CATEGORIA_TIPOLOGIA
## 1 ESSIDT
## 2 ESSIDT
## 3 ESSDIT
## 4 ESSID
## 5 ESSID
## 6 ESSID
## CATEGORIA_TIPOLOGIA
## 1 Sin atención médica en forma periódica (menor a 3 veces por semana)
## 2 Con atención médica general por lo menos 3 días de la semana
## 3 Con atención médica diaria y con especialidades y/o otras profesiones
## 4 Laboratorio de Análisis Clínicos
## 5 Laboratorio de Análisis Clínicos
## 6 Laboratorio de Anatomía patológica
## CODIGO_DEPENDENCIA DEPENDENCIA ORIGEN_FINANCIAMIENTO CODIGO_PROVINCIA
## 1 21 Provincial Público 38
## 2 21 Provincial Público 38
## 3 23 Privado Privado 26
## 4 23 Privado Privado 26
## 5 23 Privado Privado 26
## 6 23 Privado Privado 26
## PROVINCIA CODIGO_DEPARTAMENTO DEPARTAMENTO CODIGO_LOCALIDAD
## 1 Jujuy 84 Susques 38084030000
## 2 Jujuy 63 San Pedro 38063090000
## 3 Chubut 21 Escalante 26021030018
## 4 Chubut 35 Futaleufú 26035030000
## 5 Chubut 35 Futaleufú 26035030000
## 6 Chubut 35 Futaleufú 26035030000
## LOCALIDAD CP DOMICILIO TE1
## 1 EL TORO 4411 AVENIDA LAVALLE ESQUINA 24 DE OCTUBRE
## 2 LA MENDIETA 4522 CARLOS GARDEL S/N° - LOTE BARRO NEGRO
## 3 COMODORO RIVADAVIA 3240 Avenida Rivadavia 633 0297-4475310
## 4 ESQUEL 9200 25 de Mayo 498 2945452050
## 5 ESQUEL 9200 Ameghino 1033 02945-452155
## 6 ESQUEL 9200 Belgrano 347
## TE2 TE3 TE4 ESTADO_GEO geometry
## 1 SINCONFIRMAR POINT (-66.82709 -23.20853)
## 2 SINCONFIRMAR POINT (-64.93175 -24.30814)
## 3 CONFIRMADO POINT (-67.48355 -45.86213)
## 4 294515503809 CONFIRMADO POINT (-71.31971 -42.91505)
## 5 CONFIRMADO POINT (-71.31803 -42.91577)
## 6 SINCONFIRMAR POINT (-71.31992 -42.91788)Como ven el objeto salud pasó de ser un dataframe a un sf dataframe, o sea un dataframe geográfico. Y también podemos ver todas las características del objeto. El tipo de geometría es punto, tiene 2 dimensiones y el CRS es 4326 el mismo que el de radios.
Ahora podemos mapear. Y lo hacemos igual que con el dataset de barrios y radios.
ggplot()+
geom_sf(data=salud)
No tenemos los límites de Argentina pero claramente los puntos de los efectores parecen distribuirse en todo el país. Pero nosotros queremos trabajar solo con datos de la Ciudad de Buenos Aires! Una opción es filtrar a través de un atributo “no-geográfico” como por ejemplo el nombre de la Provincia, pero en este caso no tenemos esa información!
A no desesperarse! Una de las bondades de los datos geográficos es que podemos hacer consultas espaciales. En este caso vamos a pedirle a R que espacialmente nos filtre y seleccione todos los efectores de salud que caen dentro de algún barrio de Buenos Aires. Para eso usamos st_intersection una de las funciones espaciales de la librería sf.
Cargamos un nuevo dataframe con el polígono de la Ciudad de Buenos Aires.
caba <- read_sf("https://raw.githubusercontent.com/martoalalu/clase-geo-salud/clase-2021/data/caba.geojson")
ggplot(caba) +
geom_sf()
Queremos quedarnos solo con los efectores de salud (puntos) que caen dentro de la Ciudad de Buenos Aires (polígono). Entonces indicamos nuestro objeto geográfico a filtrar (efectores de salud) y luego el dataframe con el que filtrar (poligono de Buenos Aires).
#Nos quedamos solo con los efectores de salud que están dentro de barrios
salud_caba <- st_intersection(salud, caba)La consulta espacial dio como resultado un nuevo dataframe que nos indica que en la Ciudad de Buenos Aires hay 1165 efectores de salud.
Ahora sí, volvamos al mapa. Al igual que en cualquier gráfico de ggplot() podemos combinar capas. Vamos a agregar una capa el polígono de la Ciudad de Buenos Aires detrás para tener la referencia.
ggplot()+
geom_sf(data=caba)+
geom_sf(data=salud_caba)
A priori la distribución es parecida a la de densidad poblacional!
Ejercicio. Supongamos que querramos hacer un mapa que muestre la capacidad de respuesta del sistema de salud porteño, o sea sólo aquellos efectores que tengan internación. ¿Cómo lo hacemos? Va tip, el campo TIPOLOGÍA nos va a ayudar.
Volvamos al mapa, pero profundizando el análisis. Vamos a pintar cada efector según la columna ORIGEN_FINANCIAMIENTO el cual indica si es público o privado para ver cómo es la distribución espacial.
ggplot()+
geom_sf(data=caba)+
geom_sf(data=salud_caba, aes(color=ORIGEN_FINANCIAMIENTO))
Para emprolijar podemos hacer facetar el mapa, al igual que con cualquier visualización de ggplot().
ggplot()+
geom_sf(data=caba)+
geom_sf(data=salud_caba, aes(color=ORIGEN_FINANCIAMIENTO))+
facet_wrap(~ORIGEN_FINANCIAMIENTO)
4 Más allá del choropleth
Hagamos otros tipos de mapas. Tener ubicaciones exactas, es decir objetos que son puntos, nos permiten hacer otros tipos de mapas como de calor (heatmap), de burbujas (bubblemap) o de densidad (densitymap), entre otros.
Dejemos de usar un rato la librería sf y vamos a considerar nuestro dataframe de efectores de salud en Buenos Aires como cualquier otro no geográfico.
Antes con st_as_coord() habíamos transformado las columnas de latitud y longitud en una geometry ahora vamos a reestablecer esas columnas, o sea vamos a hacer que dejen de ser columnas geográficas y ser sólo “numéricas”.
#Reestablecemos las columnas lat y long con sus coordenadas
salud_caba$long <- st_coordinates(salud_caba)[,1]
salud_caba$lat <- st_coordinates(salud_caba)[,2]
#Le decimos a R que salud_caba_intern deja de ser un objeto geográfico
st_geometry(salud_caba) <- NULL
class(salud_caba)## [1] "data.frame"Ahora salud_caba es un dataframe normal que tiene una columna latitud y otra longitud que expresan la posición de los elementos en la superficie terrestre. O sea que si hacemos un
Volvamos a ggplot() y hagamos un geom_point() tradicional.
ggplot() +
geom_sf(data=caba) +
geom_point(data = salud_caba, aes(x = long, y = lat,color = ORIGEN_FINANCIAMIENTO))
Ven que queda igual que si fuera un objeto geográfico! Aprovechemos nuevas funciones de ggplot() para hacer nuevos gráficos que “simulen” ser un mapa. Un mapa de densidad de puntos que muestre las zonas con mayor concentración de efectores.
ggplot() +
geom_sf(data=caba) +
geom_bin2d(data = salud_caba, aes(x = long, y = lat), bins = 50) +
scale_fill_viridis_c()+
labs(title = "Concentración de efectores de salud",
fill = "Cantidad")
Si no nos gusta que sean cuadrados pueden ser hexagonos!
ggplot() +
geom_sf(data=caba) +
geom_hex(data = salud_caba, aes(x = long, y = lat), bins = 50) +
scale_fill_viridis_c()+
labs(title = "Concentración de efectores de salud",
fill = "Cantidad")
O bien, las zonas de calor pueden ser más estilo “curvas de nivel”.
ggplot() +
geom_sf(data=caba) +
stat_density2d(data = salud_caba, aes(x = long, y = lat, fill=stat(level)),geom="polygon")+
scale_fill_viridis_c()
5 Calculando distancias
Calcular distancias es una de las funciones espaciales más útiles ya que permite ver con claridad cómo la accesibilidad a distintos bienes y servicios se distribuye desigualmente en el espacio. En tiempos de coronavirus y cuarentena ésta cuestión cobró mayor relevancia ya que al limitarse la movilidad de las personas se cristalizan los patrones espaciales, haciendo que una parte de la población tenga un acceso privilegiado a bienes y servicios mientras que otra queda se ve obligada a trasladarse más para conseguir lo mismo (y por ende se expone más al virus).
Una de las cuestiones que surgió fue la limitación para moverse a no más de 500 metros a la redonda. Analicemos esta medida en términos de acceso a un espacio verde, algo que la Organización Mundial de la Salud considera vital para el bienestar de la población. ¿Qué porción de la población de la Ciudad de Buenos Aires está a 500 metros de una plaza?
Hagamos un choropleth map que muestre con mayor intensidad aquellos radios censales que están más cerca de una plaza.
Para esto usaremos 2 datasets: radios censales y ubicación de espacios verdes, del portal de datos abiertos del Gobierno de la Ciudad de Buenos Aires.
5.1 Explorando los datos
#Cargamos los datos
radios <- st_read('https://github.com/martoalalu/clase-geo-salud/raw/master/data/caba_radios.geojson')## Reading layer `caba_radios' from data source `https://github.com/martoalalu/clase-geo-salud/raw/master/data/caba_radios.geojson' using driver `GeoJSON'
## Simple feature collection with 3554 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -58.53092 ymin: -34.70574 xmax: -58.33455 ymax: -34.528
## Geodetic CRS: WGS 84plazas <- st_read("https://github.com/martoalalu/clase-geo-salud/raw/master/data/espacios-verdes-catastrales.geojson")## Reading layer `espacios-verdes-catastrales' from data source `https://github.com/martoalalu/clase-geo-salud/raw/master/data/espacios-verdes-catastrales.geojson' using driver `GeoJSON'
## Simple feature collection with 1402 features and 26 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -58.53175 ymin: -34.70481 xmax: -58.33983 ymax: -34.52654
## Geodetic CRS: WGS 84ggplot() +
geom_sf(data=plazas, fill="green")
Tenemos 1402 espacios verdes. Veamos qué categorías hay en el campo TIPO_EV.
unique(plazas$TIPO_EV)## [1] "PLAZA" "PLAZOLETA"
## [3] "JARDÃ\u008dN" "PARQUE"
## [5] "JARDÃ\u008dN BOTÃ\u0081NICO" "PASEO"
## [7] "CANTERO CENTRAL" "CANTEROS CENTRALES"
## [9] NA "GENERAL PAZ"
## [11] "ROTONDA" "DELLEPIANE"
## [13] "PATIO RECREATIVO Nº12" "PATIO PORTEÑO"
## [15] "PATIO RECREATIVO" "CANTEROS"
## [17] "CANTERO" "JARDINES"
## [19] "BOSQUE" "ESPACIO PÚBLICO"
## [21] "PLAZOLETAS" "RESERVA ECOLÓGICA"
## [23] "PATIO DE JUEGOS" "JARDÃ\u008dN ZOOLÓGICO"
## [25] "PARQUE DEPORTIVO" "ESPACIO VERDE"Vemos que hay espacios verdes que son canteros, plazoletas o en autopistas. Nuestro objetivo es ver el acceso a espacios verdes en los que las personas puedan distenderse y realizar algún tipo de actividad, por eso vamos a quedarnos sólo con las plazas, parques, jardín botánico, bosque, reserva ecológica, parque deportivo y espacios verdes.
plazas <- filter(plazas,TIPO_EV %in% c("PLAZA","PARQUE","JARDÍN BOTÁNICO","BOSQUE","RESERVA ECOLÓGICA","ESPACIO VERDE","PARQUE DEPORTIVO", "ESPACIO PÚBLICO"))
ggplot() +
geom_sf(data=plazas, fill="green")
5.2 Del polígono al centroide
Para calcular distancias necesitamos hacerlo de un punto exacto a otro. Necesitamos transformar nuestros polígonos (plazas y radios) en puntos exactos, para eso vamos a crear 2 nuevas capas que tengan el centroide de cada plaza y radio. Con st_point_on_surface nos aseguramos que el punto caiga adentro del polígono.
plazas_c <- st_point_on_surface(plazas)
radios_c <- st_point_on_surface(radios)Volvamos al mapa.
ggplot() +
geom_sf(data=radios_c, size=0.5) +
geom_sf(data=radios, fill=NA) +
geom_sf(data=plazas, fill=NA) +
geom_sf(data=plazas_c, color="lightgreen")
5.3 De un centroide a otro
Perfecto. Ahora bien, lo que queremos hacer es para cada centroide del radio censal encontrar el centroide de la plaza más cercana, calcular su distancia y agregarla como columna al dataframe de radios.
Se hace en pocas líneas pero vamos a explicarlo paso a paso
Primero tenemos que pedirle a R que por cada centroide de radio censal busque el centroide de la plaza más cercana. Esto lo hacemos utilizando la función st_nearest_feature.
st_nearest_feature(radios_c,plazas_c)[1:5]## [1] 283 12 2 153 333radios_c[1,]## Simple feature collection with 1 feature and 8 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -58.36622 ymin: -34.58877 xmax: -58.36622 ymax: -34.58877
## Geodetic CRS: WGS 84
## RADIO_ID BARRIO COMUNA POBLACION VIVIENDAS HOGARES HOGARES_NBI AREA_KM2
## 1 1_1_1 RETIRO 1 336 82 65 19 1.798997
## geometry
## 1 POINT (-58.36622 -34.58877)plazas_c[283,]## Simple feature collection with 1 feature and 26 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -58.37197 ymin: -34.59004 xmax: -58.37197 ymax: -34.59004
## Geodetic CRS: WGS 84
## OBJECTID SECCION MANZANA PARCELA SMP TIPO_EV NOMBRE_EV
## 283 17 003 068 000 003-068-000 PLAZA CANADÃ\u0081
## UBICACION
## 283 ANTARTIDA ARGENTINA, AV. - SAN MARTIN - MARTINEZ ZUVIRIA, GUSTAVO DR. - RAMOS MEJIA, JOSE M., DR., AV.
## OBS BARRIO COMUNA SUPERFICIE LEY FECHA_LEY ORDENANZA
## 283 DENOMINACIÓN VERIFICADA RETIRO 1 15123 0 <NA> 16976
## FECHA_ORD DECRETO FECHA_DEC BOLETIN_OF FECHA_BO FUENTE1 FUENTE2
## 283 30/12/1960 22257 23/12/1960 11505 1961-01-09 <NA> PLANO INDICE
## FUENTE3 FUENTE4 DOMINIO NIVEL geometry
## 283 <NA> <NA> GCBA 0 POINT (-58.37197 -34.59004)Lo que nos devuelve es una lista en la que asocia el index de cada elemento de radios_c con el index más cercano de plazas_c. Así sabemos que el primer radio censal (index = 1, Radio_id = 1_1_1) tiene como plaza más cercana a la ubicada en el index = 283 (Plaza Canadá).
Lo que vamos a hacer es calcular distancia elemento por elemento y para ello necesitamos tener “parejas” de radios y plazas, que por cada radio haya sí o sí una plaza asociada. Esto es sumamente útil ya que R va a agarrar el primer elemento de radios_c y va a calcular la distancia a su plaza asociada por orden, es decir uno a uno.
Entonces el siguiente paso es tener un dataframe que no sólo tenga el index de la plaza más cercana a cada radio sino tener la información geográfica (especialmente la columna geometry) para poder luego calcular la distancia propiamente dicha.
Vamos a traer la información de las plazas asociada al resultado que tuvimos con st_nearest_feature().
plazas_c[st_nearest_feature(radios_c, plazas_c),][1:5]## Simple feature collection with 3554 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -58.52777 ymin: -34.69755 xmax: -58.35335 ymax: -34.54134
## Geodetic CRS: WGS 84
## First 10 features:
## OBJECTID SECCION MANZANA PARCELA SMP
## 283 17 003 068 000 003-068-000
## 12 49 005 022 024a 005-022-024a
## 2 5 001 001 018a, 020b 001-001-018a, 020b
## 153 533 001 059A 000A,000B 001-059A-000A,000B
## 333 42 005 075A,075B,075C 000 005-075A,075B,075C-000
## 333.1 42 005 075A,075B,075C 000 005-075A,075B,075C-000
## 333.2 42 005 075A,075B,075C 000 005-075A,075B,075C-000
## 12.1 49 005 022 024a 005-022-024a
## 333.3 42 005 075A,075B,075C 000 005-075A,075B,075C-000
## 206 1045 012 008C 000 012-008C-000
## geometry
## 283 POINT (-58.37197 -34.59004)
## 12 POINT (-58.38813 -34.60519)
## 2 POINT (-58.37828 -34.60782)
## 153 POINT (-58.36943 -34.60478)
## 333 POINT (-58.38162 -34.60371)
## 333.1 POINT (-58.38162 -34.60371)
## 333.2 POINT (-58.38162 -34.60371)
## 12.1 POINT (-58.38813 -34.60519)
## 333.3 POINT (-58.38162 -34.60371)
## 206 POINT (-58.38767 -34.60916)nrow(plazas_c[st_nearest_feature(radios_c, plazas_c),])## [1] 3554Presten atención a la longitud de filas que devuelve este filtro, es 3554, la misma cantidad de observaciones que radios_c. Si bien estamos usando el dataframe plazas_c lo que le pedimos a R es que traiga los valores asociados sobre la función st_nearest_feature, que devolvía para cada index de radios_c el index de plaza_c más cercana. Asi el index = 1 del df de radios tiene como plaza más cercana al primer valor del calculo reciente (el index 283 de plazas, el object_id = 17, la plaza Canadá).
Perfecto, ya tenemos un nuevo dataframe de plazas cercanas a los radios el cual está ordenado de modo tal que el primer elemento de radios_c tiene como plaza más cercana al resultado de este filtro.
Lo que nos falta es saber la distancia en sí, lo cual hacemos con st_distance. El primer elemento de la función es el punto de origen, en este caso el centroide del radio censal, y el segundo elemento es el centroide de la plaza más cercana, el cual obtuvimos en los pasos anteriores. Luego le decimos a R que aplique la función elemento por elemento, es decir que calcule la distancia del primer elemento de radios_c al primer elemento de plazas_c[st_nearest_feature(radios_c, plazas_c),], que como sabemos es su plaza más cercana.
Veamos lo que devuelve esta función para 2 registros nada más.
st_distance(radios_c, plazas_c[st_nearest_feature(radios_c, plazas_c),], by_element = TRUE)[1:2]## Units: [m]
## [1] 545.6843 270.0224Perfecto! Ya tenemos el listado con las distancias calculadas desde el centroide de los radios al centroide de las plazas! Recordemos que el objetivo es hacer un choropleth map con la distancia de los radios a las plazas, y para eso necesitamos polígonos, no puntos. Necesitamos que los centroides de los radios vuelvan a ser polígonos.
Aca vamos a hacer un truquito. Como el largo del df de poligonos de radios es el mismo que el de los centroides (porque fue una transformación del primero) vamos a calcular las distancias entre el centroide de los radios y el de las plazas pero el resultado de esta operación lo vamos a agregar como columna en el dataframe de poligonos de radios censales. Como tiene el mismo largo R ni se va a enterar y nos va a dejar agregarlo sin problema!
(Acá tengan un toque de paciencia!)
radios <- radios %>% #Tomamos de base el df poligonos de radios censales
mutate(distancia=st_distance(radios_c, plazas_c[st_nearest_feature(radios_c, plazas_c),], by_element = TRUE)) %>% #Agregamos columna nueva con los valores del calculo de distancia entre centroidoes ;)
mutate(distancia=as.numeric(distancia)) #Convertimos a númeroMiremos como quedó el dataframe de radios.
head(radios)## Simple feature collection with 6 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -58.38593 ymin: -34.60846 xmax: -58.35054 ymax: -34.57865
## Geodetic CRS: WGS 84
## RADIO_ID BARRIO COMUNA POBLACION VIVIENDAS HOGARES HOGARES_NBI
## 1 1_1_1 RETIRO 1 336 82 65 19
## 2 1_12_1 SAN NICOLAS 1 341 365 116 25
## 3 1_12_10 SAN NICOLAS 1 296 629 101 1
## 4 1_12_11 SAN NICOLAS 1 528 375 136 7
## 5 1_12_2 SAN NICOLAS 1 229 445 129 16
## 6 1_12_3 SAN NICOLAS 1 723 744 314 104
## AREA_KM2 geometry distancia
## 1 1.79899705 MULTIPOLYGON (((-58.37189 -... 545.6843
## 2 0.01856469 MULTIPOLYGON (((-58.38593 -... 270.0224
## 3 0.04438025 MULTIPOLYGON (((-58.37879 -... 259.1591
## 4 0.36634000 MULTIPOLYGON (((-58.36733 -... 239.2970
## 5 0.01836301 MULTIPOLYGON (((-58.38454 -... 242.6532
## 6 0.03672540 MULTIPOLYGON (((-58.38154 -... 186.1420Excelente, tenemos una columna nueva (distancia) con un número que nos va a permitir hacer nuestro querido choropleth.
Al mapa directo entonces.
ggplot() +
geom_sf(data = radios, aes(fill = distancia), color = NA) +
scale_fill_viridis_c() +
labs(title = "Distancia a plaza más cercana",
subtitle = "Ciudad Autónoma de Buenos Aires",
fill = "Distancia a plaza más cercana")+
theme_void()
Genial! Si queremos podemos quedarnos sólo con los radios censales que están a menos de 500 metros.
ggplot() +
geom_sf(data = filter(radios,distancia < 501), aes(fill = distancia), color = NA) +
scale_fill_viridis_c() +
labs(title = "Radios censales a menos de 500 metros de una plaza",
subtitle = "Ciudad Autónoma de Buenos Aires",
fill = "Distancia a plaza más cercana")+
theme_void()
Ok. Muy lindo el mapa pero ¿de cuántas personas estamos hablando?
radios %>%
st_drop_geometry() %>%
group_by(distancia > 500) %>%
summarise(total = sum(POBLACION),
pct= total / sum(radios_c$POBLACION)) ## # A tibble: 2 x 3
## `distancia > 500` total pct
## <lgl> <dbl> <dbl>
## 1 FALSE 1849086 0.640
## 2 TRUE 1041065 0.360El 64% de la población de Buenos Aires está a menos de 500 metros de una plaza, mientras que el 36% no.
Por último podemos hacer un histograma para ver cómo se distribuye la distancia.
ggplot(radios) +
geom_histogram(aes(distancia))
Y lo podemos facetar por comuna…
ggplot(radios) +
geom_histogram(aes(distancia)) +
facet_wrap(~COMUNA, scales="free_y")
6 Repaso
Y asi llegamos al final de la segunda clase donde tuvimos un pantallazo de las operaciones espaciales más comúnes:
- Área -> st_area
- Perímetro -> st_length
- Intersección -> st_intersection
- Convertir polígonos en puntos -> st_point_on_surface
- Encontrar el objeto más cercano -> st_nearest_feature
- Distancia -> st_distance